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Confexf. Measurements of cosmic ray fluxes by the PAMELA and CREAM experiments show unexpected spectral features between 
200 GeV and 100 TeV. They could be due to the presence of nearby and young cosmic ray sources. This can be studied in the myriad 
model, in which cosmic rays diffuse from point-like instantaneous sources located randomly throughout the Galaxy. 
Aims. To test this hypothesis, one must compute the flux due to a catalog of local sources, but also the error bars associated to this 
quantity. This turns out not to be as straightforward as it seems, as the standard deviation is infinite when computed for the most 
general statistical ensemble. The goals of this paper are to provide a method to associate error bars to the flux measurements which 
has a clear statistical meaning, and to explore the relation between the myriad model and the more usual source model based on a 
continuous distribution. 

Methods. To this end, we show that the quantiles of the flux distribution are well-defined, even though the standard deviation is 
infinite. They can be used to compute 68 % confidence levels, for instance. We also use the fact that local sources have known 
positions and ages to reduce the statistical ensemble from which random sources are drawn in the myriad model. 
Results. We present a method to evaluate meaningful error bars for the flux obtained in the myriad model. In this context, we also 
discuss the status of the spectral features observed in the proton flux by CREAM and PAMELA. 
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1. Introduction 

In the 1 GeV to 100 TeV energy range, cosmic ray (CR) nu- 
clei that reach the Earth have a Galactic origin. They were ac- 
celerated in sources, the nature of which is still a subject of 
discussion and research, supernova driven shock waves being 
a good candidate. They subsequently reach the Earth, after dif- 
J> fusing in the Galactic magnetic field. The exact locations and 
Q\ ages of the supernova explosions are not known, and properties 
OO of Galactic cosmic rays are often studied under the assumption 
that sources form a distribution which is continuous, in space as 
in time. In the standard model of CR propagation, sources are 
■^j- modelled as a jelly which extends inside the Galactic disk and 
steadily injects particles in the interstellar medium (ISM). This 
picture has proven to be quite successfull so far. The CR spectra 
t-H are actually dominated by sources which are distant from Earth 
^ (for instance see Fig. 8 in ( jTaillet et al.||2004") >) and the contin- 
• i-h uous hypothesis is expected to be valid. On the other hand, this 
approach may fail, when one considers the high-energy end of 
?-j the CR spectra, mainly produced by young sources. The flux at 
^ high energy could well be dominated by a handful of individual 
sources, the positions and ages of which have a strong influence 
on the spectrum. This is of tremendous importance as the latest 
measurements, such as those of CREAM ( |Yoon et "aT7||201 1] > or 
PAMELA ( |Adriani et al.|201l| l, exhibit spectral features which 
are difficult to understand in a continuous and steady-state mod- 
elling of CR sources. 

In this article, we study the relation between the mean- 
field approach and the so-called "myriad model" ((Higdon et al. 
|2003) l) where CR primary species, such as protons or helium 
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nuclei, diffuse from point-like sources distributed randomly in 
space and time. In the myriad model, the flux of cosmic rays 
is a random variable. When averaged over many realisations of 
the possible positions and ages of the sources, the myriad model 
gives the same results as the continuous model, as expected. 

However, we expect fluctuations from the average flux, from 
one realization to the other. If these fluctuations are small, then 
the myriad model predictions of the flux are robust and similar 
to that of the continuous model. In that case, the spectral features 
seen by CREAM and PAMELA should be explained within the 
continuous model. Invoking a distorted spectral injection spec- 
trum at the sources (see for instance fYuan et aL] ( |201 l| l) or a non- 
standard behaviour of the Galactic diffusion coefficient with en- 
ergy are two possible, although unlikely, solutions. On the other 
hand, if the statistical fluctuations of the myriad model are large, 
it could be easier, and to some extent more natural, to explain 
the CREAM and PAMELA anomalies in terms of the particu- 
lar history and morphology of the population of sources within 
which we live. If so, there is little reason to trust anymore the 
continuous model, on which most of CR studies are neverthe- 
less based. The "real" sources could generate CR fluxes at the 
Earth quite different from those predicted in the conventional 
approach. Unfortunately, we do not know where and when the 
"real" supernova explosions took place. We are reduced to rely 
on statistics in order to gauge how probable are CR fluxes far 
from the values predicted by the continuous model. 

The amplitudes of the statistical fluctuations of the myriad 
model have actually been recently studied through the standard 
deviation of the random variable associated to the CR proton flux 
(see e.g. (Blas i & Amato|20"TT i). This quantity turns out te be in- 
finite, when computed without precautions. This unexpected re- 
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suit is, a priori, a severe blow against the mean-field approach on 
which most of the public codes of CR Galactic transport, such 
as galproiQ DRAGOt^ or usin^] are based. These codes would 
have to be entirely modified. Abandoning the conventional CR 
model would have also drastic consequences on the problem of 
the astronomical dark matter (DM). The nature of this compo- 
nent, which contributes substantially to the mass budget of the 
universe, is still an unresolved issue. Weakly interacting mas- 
sive particles (WIMP) have been invoked as a plausible solution. 
These putative DM species are expected to annihilate inside the 
Galactic DM halo, producing in particular rare antimatter cosmic 
rays such as positrons and antiprotons. Distortions in the spectra 
of these particles would be an indirect probe of the presence of 
WIMPs inside the Milky Way. Theoretical predictions have been 
based so far on the continuous CR model. They need to be en- 
tirely revisited should the myriad model replace the conventional 
scheme. 

This article is devoted to the myriad model and to the is- 
sue arising from the infinite variance of the primary CR flux. 
We suggest two ways to solve this problem. First, we show that 
the statistical distribution of the flux is such that the confidence 
intervals are finite and well-defined, even though the variance 
is infinite. We study how the width of these intervals depends 
on the CR propagation model. Second, we notice that the diver- 
gence of the variance is due to the sources that are very close 
to the Earth and very young. In our Galaxy, theses objects are 
known actually, and the region in which we have catalogs of the 
sources should be excluded from the statistical ensemble of the 
myriad model, and rather treated as known sources. This paper 
explores the consequences of this approach. 



2. Propagator for diffusive propagation from 
discrete sources 

After their acceleration by supernova driven shock waves, pri- 
mary cosmic ray (CR) nuclei subsequently propagate through 
the Galactic magnetic fields and bounce off their irregularities 
- the Alfven waves. The resulting particle transport is well de- 
scribed by space diffusion with a coefficient 



D(E)=D /3, 



(1) 



which increases as a power law with the rigidity 3£ = p/Z of the 
particle. In addition, since these scattering centers move with a 
velocity V a ~ 20 to 100 km s _1 , a second order Fermi mechanism 
is responsible for some diffusive re-acceleration which turns out 
to be mostly relevant at low energy, below a few GeV. Since we 
are interested here in the excess in CR protons and helium nuclei 
measured by CREAM, and more recently by PAMELA, we will 
disregard diffusive re-acceleration as well as energy losses which 
are also negligible at high energy. In addition to space diffusion, 
Galactic convection wipes cosmic rays away from the Galactic 
disk with a velocity Vc ~ 5 to 15 km s . The master equation for 
the space and energy number density if/ = dn p /dT p of, say, CR 
protons with kinetic energy T p takes into account space diffusion 
and Galactic convection 



dif/ 
~dt 



+ d z (V c iff) — D Aif/ — q ef f = g acc - q c , 



(2) 



1 http://galprop.stanford.edu/ 

2 http://www.desy.de/ maccione/DRAGON/ 
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This equation applies to any nuclear species - protons and he- 
lium nuclei in particular - as long as the effective rate of produc- 
tion q e ff is properly accounted for. This general framework, sum- 
marized in Eq. is generally implemented within an axisym- 
metric two-zone model which is extensively discussed in Maurin 
|et al.| ( |200T| ) or |Donato et al.| ( |2001| l, and whose salient features 
are briefly recalled here. 

The CR diffusive halo (DH) is pictured as a thick wheel 
which matches the circular structure of the Milky Way. The 
Galactic disk of stars and gas - where primary CR protons and 
helium nuclei are accelerated - lies in the middle. It extends radi- 
ally 20 kpc from the center and has a half-thickness h of 100 pc. 
Confinement layers, where cosmic rays are trapped by diffusion, 
lie above and beneath that disk. The interGalactic medium starts 
at the vertical boundaries z - +L as well as beyond a radius of 
r — Rgai — 20 kpc. Notice that the half-thickness L of the DH 
is not known and possible values range from 1 to 15 kpc. The 
diffusion coefficient D is taken to be the same everywhere while 
the convective velocity is exclusively vertical, with component 



Vc(z) = V c sign(z) 



(3) 



The Galactic wind - produced by the bulk of the disk stars like 
the Sun - drifts away from them along the vertical directions, 
hence the particular form assumed here for V c . The effective 
production term q e ff takes into account the release of primary CR 
nuclei described by the positive term q acc as well as a negative 
contribution that accounts for their interactions with the inter- 
stellar gas of the disk with rate q co \. In most propagation codes, 
the injection of cosmic rays is described by a smooth function of 
space and is constant in time. In this article, we treat supernova 
explosions as point-like events. The production rate of CR nuclei 
through acceleration is given by 



(4) 



where each source i belonging to the population of super- 
novae contributes a factor qi at position X; and time f,-. The colli- 
sions of CR nuclei on the hydrogen and helium atoms of the in- 
terstellar medium (ISM) tend to deplete the high energy regions 
of their spectra. In the case of CR protons, which we will use 
throughout this section as an illustration, they act as a negative 
source term for the energy bin located at T p , with amplitude 



<7coi(x.s , t s ) = 2 h 6(z s ) T p if/(x s , t s ) 



(5) 



The ISM is distributed within an infinitely thin disk, hence the 
presence of a S(zs) function in the previous expression. The col- 
lision rate may be expressed as 



T p = V p (cTpH n H + CTpHe «He) , 



(6) 



where the densities «h and «ee have been respectively averaged 
to 0.9 and 0.1 cm 3 . The total proton-proton cross sectio n cr pii 



has been parameterized according to |Nakamura et aL 



(2010 



http://lpsc.in2p3.fr/usine 



while o-pn e is related to <x pH by thefNorbury & Townsend (2007 
scaling factor 4 2 2 ^ 3 . Similar scaling factors have been used in or- 
der to derive the CR helium nuclei collision cross sections from 
the proton case. 

2.1. Solution of the diffusion equation 

Our description of the propagation of CR protons through the 
DH, and of any other CR nucleus for that matter, relies on the 
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existence of the propagator This Green function translates the 
probability for a CR proton injected at position Xs = (xs ,ys,Zs) 
and time t$ to travel through the Galactic magnetic field until it 
reaches, at time f, an observer located at x = (x,y, z). The CR 
proton density at the Earth may be expressed as the convolution 
over space and time of the Green function & p with g acc 



(Kx, f) 



f dt s f c 

J-oo Jdh 



,(x,r <- %s,t s ) quotes, ts) 



(7) 

where q acc (^-s Js) is the CR proton injection rate at the source 
located at x s and at time t$. Notice that for a population 
of point-like sources, the injection rate q acc is given by Eq. Q 
which, once plugged in the previous expression for leads to 
the proton flux at the Earth 



4-7T 



\tffs %(%0,t = O <- X„,t = t n )Xq n \ . (8) 

I ne^ 2 j 



The Green function 8^, plays a crucial role in our analysis and is 
the ideal tool to study the effect of the discreteness of the sources. 
It is a solution of the CR transport equation 



d_% 
dt 



y+d z {v c %)-DA% + 2h 5{z) T„ % = c5 3 (x-x s ) 5{t-t s ) 



The resolution of this equation is presented in Appendix [A] 



(9) 



extend to R m!a — > oo. We obtain 



T^max^max *Jo Jlmin 



c/r 



4ttDLt 



-r 2 /4£>T ^ ( 



4'Dt 



Considering v sources per unit time and unit area in the disk, the 
average flux from N - V7i7^ ax f max sources is given by 



4nDLj Jo r 4^ 



By integrating over time (f max — > oo) and using ( Gradshteyn et al. 
( (5507) 3.471.9) 

x v - 1 e-V*-?* dx = 2 ^ ' K v (2 ViSy) , (10) 

where K v is the v-th order modified Bessel function of the second 
kind, leads to 



<Sfo> = v 



27ZTSf ste ady(>) t/r 



where Steady stands for the one source steady-state propagator, 
given by 

1 °° 

Steady (« *~ X S ) = J] ^ (p V^Td) . (1 1) 



2.2. Diffusion parameters 

The computed flux depends on the diffusion parameters. The re- 
sults will be presented for three benchmark sets of diffusion pa- 
rameters, consistent with the energy dependance of the B/C ratio 
( |Maurin et al.|200l) . They are labelled "min", "med" and "max", 
according to the value of L. The values are indicated in Table [T] 



Table 1. Diffusion parameters for the three benchmark models 



model 


AiCkpc-My- 1 ) 


S 


L (kpc) 


V c (km • s- 1 ) 


min 


0.0016 


0.85 


1 


13.5 


med 


0.0112 


0.7 


4 


12 


max 


0.0765 


0.46 


15 
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3. Mean flux 

The flux at Solar position is obtained by summing the contribu- 
tions of all the point sources. The mean value can be computed, 
once the distribution of sources is known. We first provide the 
result for a homogeneous distribution of sources, as an analytic 
expression can be found in this case. We then turn to a more 
realistic distribution. 



3.1. Homogeneous distribution of sources 

The average flux from one point source, drawn from a statisti- 
cal ensemble with distances ranging from R mm to R mdx and ages 
from fmin to f max , is given by averaging Sfi = & p over distances 
and times. To get an order of magnitude estimate, we first con- 
sider that the sources are distributed evenly in space and time in 
the Galactic disk, which is assumed to be infinitely thin and to 



Finally, using J xKo(x)dx - —xKi(x), we have 

1 In + 1 nR 

n=0 

-RmaxKi 



°° 2 ( (' 
Dj-i(2n+Y)7t \ \ 



2 L 
2n + 1 ttR„ 



(12) 



For an infinite disk, R m i n - and R mdx — > oo, so that, using 
xK\(x) — > 1 when x — > 0, 



- — y( 2 ) 2 

} ~~D l?M2n + l)n) 



n=0 

It can be shown that l/(n + 1/2) 2 = n 1 12, and finally for the 
infinite disk, 

vL 
2D 

which is what is also obtained by directly solving the diffusion 
equation in steady state. The mean value of the flux from ran- 
domly distributed point sources is equal to the steady-state flux 
obtained with a continuous source distribution. This is also true 
in more general cases, thick disk, finite radius, with wind and 
spallation. 

3.2. Realistic distribution of sources 

We now consider the case of a general distribution of sources, 
with a radial distribution in the Galactic disk, as well as a distri- 
bution across the thickness of the disk. It is quite well accepted 
that up to energies corresponding to the knee, cosmic rays ac- 
celerators are supernova remnants (SNR). Unfortunately, super- 
novae are pretty rare events and their spatial distribution is dif- 
ficult to measure accurately. However, as pulsars are created in 
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SNR and are easier to detect, it is a fair asumption that the cos- 
mic ray sources follow the pulsar distribution. In this work we 
suppose that the radial profile of sources follows the pulsar dis- 
tribution given in Yusifov & Kiiguk (2004), 

, 1.64 
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Then we can write the mean value as 
A 



<Sfo> = N 



fr(r s ) = 



r s + 0.55 
8.5+0.55 



exp {-4.01 



r s - 8.5 
8.5+0.55 



<J ft 

(17) 

This result is integrated numerically over the two coordinates r s 
and 6 S . The mean flux is then given by 

($>) = (^N)^-g(E) 
4n 

Where g(E) is the source energy injection term. We checked that 
this calculation gives the same result as the steady state model. 

In"thi7sVc^onVz,r^Tar5Vrefer to"&eToritlon"and""age of In the ne f se f io "' we P erfom the calculation of the variance 

associated to the flux. 



where distances are expressed in kpc. The distribution along the 
Z axis is given by 

tt \ / W 
f s (z s ) = exp 

where zo is set to the half-thickness of the Galactic disk, zo = h. 



sources. The solar system is located at r Q — 8.5 kpc, 8-0 and 
z — 0. Moreover we define p(r s , S ) as the distance from a source 
to the Solar System: 

p(r s , S ) = ^(r Q - r s cos S ) 2 + r] sin 2 S 
We evaluate the scatter of the flux that arises when we consider 



4. Statistical analysis : variance of the flux 

The flux obtained in the myriad model depends on the exact po- 
sitions and ages of all the sources. It is bound to be different 

from the mean values computed above, at some level. How dif- 
point-like sources with positions and ages following a given ferent from ±e mean value are likely t0 be me ty p ical fluxes in 

the myriad model? In order to address this question, the usual 
method is to compute the standard deviation associated to the 
flux, considered as a random variable. 



probability distribution. The sources are supposed to be inde- 
pendent, so the mean value of the propagator coming from all 
the N sources is just given by : 

<^v) =Nx(& l ) = NA j 2nr s dr s j d6 s j dz s j dt s 
feM Mr,) f z (z s ) f,(t s )% (r„ 6 S , t s ) 



(13) 



Assuming cylindrical symmetry and a uniform age distribution, 
up to a maximum age T max , the normalization factor A is given 
by, 

- = 2n T max I rf r (r) f z (z)drdz 

™ ^Galaxy 

We find (the antisymmetric part of the propagator vanishes at 
Z = 0) 



4.1. Realistic distribution of sources 

The variance of the propagator for one source is defined as 



cr 2 = 



,) 2 



1 - wi / - WIj 
The variance for N independent sources is given by 
o-l = Ncr\ 

We must compute the average value of ^ 2 , given by 



d 



N J r s dr s dz s dt s d9 s A f r (r s ) exp 



l£vl 
Zo 



4nDt 



exp 



-VAzs 
2D 



, P 2 (r s ,0 s ) , 

exp( -^r } 



/dr s dz^dt s d6 s Ar s f r ( r s) exp (-— 
\ zo 



s 



exp (-(a„)t s ) 
C n 



sin(k„L) sin(k n (L - \z s \)) 



x 



4nDt s 

CO CO 



exp 



Vc\z s \ 
D 



exp 



P 2 (r s ,0s) 
2Dt s 



(18) 



exp (-(a„ + a m )t s ) 



We separate the integration and define 



sm(k„L) sin(k m L) sin{k„(L - \z s \)) &in(k m (L - \z s \)) 



{J Z )n = 



sin(k n L) f L ( \z s \\ ( V c \zs\\ ..... 



\z s \)) 
(14) 

The integration is easy using the exponential form of sinus and 
leads to 



At fixed m and n, we define 
sin(A:„L) sin(A: m L) 



^7 }n,m — 



C C 



{Jz)n = 



sin(k„L) 



C„ (V c /2D + l/h) 2 +k 2 



sm(k„(L - |z s |)) sin(k m (L - \ Zs \)) 
(19) 



^ + - J sin(k„L) - kn (cos(k„L) - e - L (ffi + i)) 



(15) 



\\2D h, 

The integration over time (f max — > oo) is made using Eq.[l0|: 

p\r s ,6 s )\ 



r°° i 

Jo 4nDt s 



exp - 



1 



4Dt s 



exp(-a n t s )dt s 



2Dn 



K \p(r s ,e s ). 



(16) 



This is easily computed, and we can integrate the time distribu- 
tion to get the variance from sources of all ages : 

, )„ >m dt s dr s dO s r s f r {r s ) exp(-(ar n + a m )t s ) 

n,m 

2p 2 (r s ,6 s )\ 



1 

l6n 2 D 2 t 2 



exp 



4Dt, 



(20) 
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Using Eq.fTUland considering ? max — > oo 



exp 



2D?. S 



exp(-(o'„ + a m )f s ) 



= ^2(a„ + or m )D 



1 1 



8tt 2 D 2 p 



2(ff„ + a m ) 



D 



(21) 



we find 



n.lll 

x I d8. 



' 8/r " f> " J " 

X 



P(r s , S ) 



p{r s ,6 s )J2 



Mr,) 
(22) 



We are interested in the behaviour of this quantity at the lower 
bound in r. When r y/2(a m + a„)/D «: 1 , the property K\ (x) — > 
1 jx as x — > can be used to show that the above integral diverges 
as ln(/? m ; n ) as the lower bound R m \ n of the spatial integral goes to 
zero. The spatial distribution probability is f r ~ 1 when r — > 0, 
so that 

<^, 2 > ~ ^ 2<^ 2 )„, m ln^ mm (23) 



4.2. Divergence of the variance: how bad is it? 

The mean value, computed considering r s from r s = to r s = 20 
kpc and f s from t , = to f s = oo is equal to the flux obtained 
in the steady-state model, as commonly obtained in propagation 
theories. But as we have just seen, the calculation related to the 
standard deviation of this quantity diverges, when we allow the 
possibility of sources that are close (r s — > 0) and young (f s — > 0), 
as already noticed before (( |Blasi &~ Amato 20TTJ). 

The standard deviation is usually interpreted as the typical 
spread of the random values around the mean, and a large stan- 
dard deviation could be interpreted as if the actual value of the 
flux had a disturbingly large probability to be very far from the 
mean value. 

One could argue that the problem we considered is physi- 
cally irrelevant, as we know for sure that there is no supernova 
remnant with zero age and null distance to the Earth. One can 
impose a lower cut-off in ages and distances, based on obser- 
vations. However, even with reasonable values for the cut-off, 



the variance can be finite but large (see section 4.3 I. One can 
also adjust the cut-off in order to eliminate the very rare events 
that make the standard deviation very large, without contributing 
significantly to the the mean value. This is the approach adopted 
by |Blasi & Amato](|2011[ ). These authors use a cut-off given by 
Wi = Rmax/ ^4vD(E). This condition gives a fair order of mag- 
nitude of the spread of the values around the mean. It is difficult 
though to interpret it in rigourous statistical terms. Actually, the 
value of the variance depends quite strongly on the exact po- 
sition of the cut-off, as featured in Fig. [TJ where the evolution 
of the mean and the standard deviation of the total proton flux 
derived with the max model at 10 TeV is presented. Choosing 
tmm /2 or 2 x f m ; n rather than f m ; n , for instance, has a small effect 
on the mean but has a drastic effect on the standard deviation. 
The value chosen by Blasi & Amato ( |201 l| l is indicated by an 
arrow on the figure. The mean obtained with this cut-off is about 
10 % smaller than the true mean. A lower cut-off would give a 
more precise mean value, but a much larger variance. 

This confusing situation, where some rare events have a very 
small contribution to the mean, but give rise to a very large 



standard deviation, is not uncommon in physics (Levy flights, 
Cauchy distribution). The total flux (f> is the sum of a myriad 
of N individual contributions ipi from single sources. Assuming 
that these contributions are not correlated with each other allows 
to relate simply the means and variances of the total and single 
source fluxes through 



(cp)=Nx (ip) and crl = N X at 



N(<p 2 ) 



N 



(24) 



We show in the appendix that the probability p(<p) of measuring 
a single source flux <p is given by p(<p) cc tp~ & / 3 for objects located 
in a thick disk and p{<p) cc ^r 7 / 3 for sources located in a thin disk, 
in the limit where <p — > oo. The standard deviation cr^ of the total 
flux is related to the integral (ip 2 ) — J ip 2 p(tp) dtp, which diverges 
as <p — > oo. For such a distribution with an infinite second mo- 
ment, the central limit theorem does not hold, at least in its usual 
form. In this case, the standard deviation is not a good estimator 
of the typical spread of the values around the mean (in particular 
the probability distribution function p(<p) is no longer Gaussian), 
and there is no use trying to regularize this quantity by applying 
cut-offs. Notice that the average value (0) is still well-defined, 
as (<p) = J ip p(ifi) d<f> is convergent, and more to the point, the 
confidence intervals of the global probability distribution p{<p) 
are well-defined too, since the integral of p{ip) is well-behaved. 
The spread of the random values of the total flux around its mean 
(<p) can then be studied by computing the quantiles associated to 
the probability distribution p(<f>), rather than using the standard 
deviation which has no clear statistical relevance in this case. 

All this discussion arises from the fact that, in the myriad 
model, there is a non-vanishing probability to find sources which 
are arbitrarily close and young. Although the variance cr^ can- 
not be defined in this case, we can still infer the global distribu- 
tion function p((p) and delineate intervals inside which the total 
fux <p is mostly expected. Another possibility is to take advan- 
tage of the solid information which has been collected on these 
young and nearby sources, whose statistical properties are dis- 
concerting. We can actually remove a space-time region around 
the Earth containing the very close and very young sources of 
our Galaxy, and replace them by a catalog, built from observa- 
tions. 

In what follows, we investigate these points in further details. 
First we study how the standard deviation cr^ can be regularized 
by applying some cut-offs related to astronomical information 
about the sources. We then focus on the confidence intervals of 
p{4>) and we show that they are well-defined, with no need of a 
particular cut-off. Finally, we take into account the astronomical 
knowledge about the local sources of cosmic rays and explore 
in greater detail how introducing such a catalog improves our 
statistical analysis. 

4.3. Regularisation by cutting off the ages 

From observations of the Solar neighbourhood, we have a good 
idea of the distribution of sources which are young and close. 
Given the age t mm of the youngest local supernova remnant, one 
can compute the mean value and standard deviation of the total 
flux (p, by applying that cut-off f m ; n to the age distribution. 

For the sake of illustration, Fig. [2] features the mean value 
{(f>) and the standard deviation cr^ of the total proton flux using a 
lower cut-off of f m ; n = 100 yr in the source distribution, for the 
min set of propagation parameters. On the same figure we have 
also plotted the data points from the CREAM and PAMELA ex- 
periments. With the chosen value for the cut-off, the standard 
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deviation at high energies remains of the same order of mag- 
nitude as the flux. The relative value of the standard deviation, 
o"^ /((/>), is fairly independent of the energy. This trend can be 
shown analytically to hold for sources located in a thin disk. In 
the framework of a purely diffusive model (no wind/spallation), 
the relative dispersion can be actually approximated by 



R 



(4>) 4LV2^ 



(25) 



This ratio does not depend on the diffusion coefficient D, hence 
it does not depend on energy. It is of order or unity for t m { n ~ 
100 yr. 

For sources distributed in a disk having a finite thickness h, 
one gets the same result as before as long as h <K y/Dt m i n . In the 
opposite limit, one finds 



cc D 



1/4 



(26) 



In both cases, the relative standard deviation does not vary much 
with energy. 

The value of cr $/(</>) depends sensitively on the chosen cut- 
off 

tmm- For a very low value of f m ; n , the standard deviation is 
very large. Conversely, for larger values of f m i n , the standard de- 
viation decreases, and the average value (<f>) starts to be affected 
by the cut-off, as shown in Fig.[T] 

4.4. The variance is infinite but the confidence levels are 
finite 

The variance of the flux being infinite does not necessarily im- 
ply that the random values are typically very different from the 
mean. To illustrate this affirmation, consider the case of a unique 
point-like steady-state source located in the Galactic disk, with 
cosmic ray diffusion taking place in a boundless space. The so- 
lution of the diffusion equation is given by ip — a/r where a is a 
constant. Assuming that this source is uniformly distributed in- 
side a disk of radius R leads to the probability distribution func- 
tion 

2nrdr 2rdr 

dp =^r = lv- (27) 

We can readily infer the mean flux 

a 2rdr 



if) 



-f 

Jo 



R 2 



2a 
~R 



(28) 



and the average value of the flux squared 



<<P 2 ) 



C R a 2 2rdr a 2 (R\ 

= 1 7^ = rAi) • (29) 



where we have introduced a cut-off value e at the lower end of 
the radial distribution, to exhibit the divergence of {if 2 ). The vari- 
ance of (p goes to infinity as e — > 0. 

However, the distribution of ip (which is what we are really 
interested in) is well-behaved. From the relation between r and 
tp, we can write dr — a dip/if 2 so that 



dp(ip) 



2r dr 2a dip 



R 2 



R 2 ip^ 



(30) 



The probability that the flux is smaller than a given value <t> may 
be expressed as 



P(< <£) 



rv(R) 
J<p(r) 



dp(<p) = 1 - 



R 2 <b 2 ' 



(31) 



Table 2. Divergence of the variance 





2D 


3D 


steady-state 
time-dependent 


(if) finite 

(if 2 ) — > oo 

(tp) finite 
{ip 2 ) -> oo 


(ip) finite 
(ip 1 ) finite 
(if) finite 

(if 2 ) — > oo 



provided that O > <t>o = o/R. Introducing the Heavyside distri- 
bution leads to 



P(> <E>) = 



R 2 ® 2 



Q 



(•-*)- 



4<D 2 l 2 



(32) 



The probability that O > 10 (tp) is only 1/400, even though the 
variance is infinite. Actually, the flux is more likely to be lower 
than the mean value, whereas one might have guessed the oppo- 
site, considering the divergence of the variance. 

When N sources are considered, the mean flux value and the 
variance are both just multiplied by N. The probability distribu- 
tion Pn(4>) for the flux can be obtained by recurrence from 



Pn(4>) 



-I 



p(ip) p N -i(4> - ip)d<p 



(33) 



These are displayed in Fig. |4] The variance still diverges. In the 
high-0 region, the flux is dominated by the contribution of a sin- 
gle source and the probability distribution is given by 



dp N 
dcf> 



N 



dp 2Na 2 



dp R 2 ^ 



(34) 



For N = 100 sources, the probability that <f> > 2 (cp) is 2.5 x 10~ 3 
and P((p > 10 (0)) is vanishingly small. 

If we now consider time-dependent sources spread homoge- 
neously inside an infinite DH with pure diffusion, the variance is 
given by the integral 



47rr dr 



1 



(4jTDt) 3 ' 2 



-r~/2Dt 



(35) 



which diverges with the lower cut-off in ages as 1 / V ? min ■ For a 
3D homogeneous distribution of steady-state sources, cr^ does 
not diverge (see Table El. For the sake of illustration, Fig. [5] 
presents histograms for the flux obtained with the propagator 
discussed in Sec. 2, at several energies. 

4.5. Quantiles of the total flux 

From now on, we show the quantiles associated to the fluxes. 
These are defined as regions with a given probability to find the 
flux. In Fig. [6] we present "deciles" (10 % quantiles), as well as 
68 % confidence intervals (CI), which are more familiar, as they 
correspond to 1-cr intervals for gaussian distributions. 

In order to compute these figures, we ran a Monte Carlo sim- 
ulation over ~ 3000 populations. We also adjusted the injection 
rate so that the mean flux follows the PAMELA data below 200 
GeV : 



3.53 x 10~ 3 x(l -■ 



in cm 2 s 'sr 'GeV 



f T . -2.5 / i ,t \ -0.42 

-(r p /2.5)""\ / 1 p \ 1 1 + 1 p ] 



10 



16 



(36) 



The spread of the fluxes around the mean is not sensitive to 
this choice. In spite of its infinite variance, the flux has a well- 
defined probability distribution as featured in Fig. [6] The possi- 
ble flux lies on a band whose thickness depends sensitively on 
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the propagation model. In particular, the width increases when 
the thickness L of the diffusive halo decreases. 

At high energy, the CREAM and PAMELA data cannot be 
explained in the med and max cases. For the min case, taking 
v = 1 century 1 would enlarge the uncertainty band, so that the 
explanation in terms of sources becomes more probable. 

5. Known local sources - The catalog 

5.1. Regulahzation using a catalog 

The situation described above occurs as long as we know noth- 
ing of the positions and ages of the CR sources. The young and 
nearby objects are responsible for the divergence of the flux vari- 
ance and potentially lead to the problems encountered above in 
the statistical analysis. However, we actually have data concern- 
ing the distribution of nearby sources, for which catalogs are 
available. A natural way then to regularize the variance is to sep- 
arate the sources in two lots. The first set contains the young 
and local sources, which can be extracted from the catalogs. The 
second group, on which little information is known, comprises 
the old or distant sources and will be treated according to the 
statistical analysis of Sec. [4 This procedure allows to regularize 
the variance of the flux in the most natural way while reducing 
its uncertainties. Following Delah aye et al.| ( [20TQ| >, we have used 
two catalogs. 

(i) The Green survey ( Green||2009| ) compiles various informa- 
tions on supernova remnants, but fails to systematically pro- 
vide their ages or the precision with which their distances from 
the Earth have been determined. A quite thorough bibliographic 
work has been summarised in the appendix of Delahaye et al. 
(2010), which we have borrowed as a complement of the Green 
catalog. In total, we have collected 27 local SNR with their ages, 
distances from the Sun and, when possible, the corresponding 
observational uncertainties. 

(ii) Pulsars are not expected to be sources of primary CR nu- 
clei. As residues of supernova explosions, they are nevertheless 
a good tracer of old SNR too old to be detected directly in radio 
waves. Moreover, being point-like objects, their distances from 
the Sun is way easier to measure. Their ages can also be esti- 
mated precisely through spin-down. After removing millisecond 
pulsars from the ATNF catalog (Mancheste r et al.|2005[ l (by se- 
lecting P > 5 x 10~ 18 ) and objects associated with known SNR, 
we are left with 157 objects with their ages and distances. 

It is fair to ask whether these 27+157 sources are representa- 
tive of the local environment. As featured in Fig.|7J the number 
of objects found in the catalogs is in good agreement with what 
can be inferred from various Galactic distributions found in the 
literature, provided that the supernova explosion rate is, on av- 
erage, approximately equal to 3 per century. At least, this is true 
within the 2 kpc nearby the Sun, for sources younger than 30,000 
years. At later ages, SNR become too dim to be detected and the 
Green survey cannot be trusted anymore. However, the catalogs 
are in rather good agreement with theoretical expectations and 
do not suffer from major biases, at least not more than the theo- 
retical models. In this work, we define the "local" region as the 
domain extending 2 kpc around the Sun, with sources younger 
than 50,000 years. A rate of 1 supernova explosion per century in 
our Galaxy is also quoted in the literature (Delahaye et al. 2010). 
If so, we would be in a locally high-density region of sources. We 
have plotted the flux considering this low rate, but will discuss 
the validity of this assumption in a forthcoming letter. 

The CR proton flux produced by the SNR and pulsars of 
our catalog is presented in Fig. [8] for the "min" CR propaga- 



tion benchmark model. In the PAMELA energy region, the flux 
is dominated by pulsars whereas SNR come into play at the en- 
ergies of CREAM. 

5.2. Is the Catalog probable? 

In order to check the plausibility of our catalog, we compare 
in Fig. [9] the flux </> cat yielded by the sources which it contains 
(blue curve) to the flux <f>i oc from a set of populations drawn ran- 
domly inside the local region (black curve). As explained above, 
we did not compute the variance of these random populations, 
but derived through Monte-Carlo simulation the confidence in- 
tervals. Assuming an explosion rate of 3 events per century in 
the Galaxy, we find that the flux </> cat yielded by the objects of 
our catalog lies within the 1-cr confidence interval surround- 
ing the mean value 0i oo . This is true for the min, med and max 
benchmark sets of CR propagation parameters. We conclude that 
the theoretical source distribution and explosion rate, which we 
have chosen here, are in rather good agreement with local real- 
istic sources. They are not necessarily representative of the en- 
tire Galaxy though. How potential differences would impact the 
relative importance of <f> cat with respect to the total flux will be 
detailed in a forthcoming letter. 

5.3. Systematic errors from the catalog 

We have also plotted the uncertainty on cat that arises from 
the ages and positions of the SNR of our catalog. These have 
been varied within the ranges allowed by observations. We did 
not consider the uncertainties on the ages and positions of pul- 
sars. The former are well determined. The rotation periods P and 
drifts in time P of pulsars are actually measured with good ac- 
curacy. The uncertainty on </> cat corresponds to the grey bands of 
Fig. [9] and Fig. 1 1 The shape of these curves is explained by 
the fact that when L is lower ("min" model), the effective range 
of diffusion in the plane of the disk is reduced and the contribu- 
tion of remote sources is smaller. The relative contribution of the 
catalog is then larger for smaller values of L. 



6. A mixed approach 

The cosmic ray flux from the Galactic sources and the associated 
uncertainties can be computed as 

0tot(£) = 0cat(#) + 0ext(£) 

where cat is the flux coming from the sources belonging to 
the catalog (closer than 2 kpc and younger than 3 x 10 5 yr) 
and (pext(EJ) is the flux coming from sources which are more 
distant than 2 kpc or older than 3 x 10 5 yr. The uncertain- 
ties on <p ca t(E) are obtained by considering the observational er- 
rors on the source parameters (see Fig. 11 1. The uncertainty on 
<pext(E) can be evaluated by computing the confidence intervals, 
as above. In some conditions, this is equivalent to computing the 
usual the standard deviation of the flux, as the variance of the 
flux located outside of the region covered by the catalog is finite. 
This is shown in Fig. 10 



7. Conclusions 

In the conventional model of Galactic CR propagation, the 
sources of primary nuclei are described as a jelly which spans 
the disk of the Milky Way and continuously injects particles in- 
side the ISM. The actual distribution is lacunary and consists 
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of point-like objects which release cosmic rays in a very short 
time. Because measurements have become very accurate, tak- 
ing into account the discreteness of the sources in our descrip- 
tion of CR propagation has become timely. In particular, the 
PAMELA ( |Adriani et all20TT) and CREAM ffoon et al.|20TT) 
observations point towards an excess in the CR proton and he- 
lium fluxes with respect to the pure power-laws predicted with 
the continuous model. Very few analyses have been devoted so 
far to the myriad model. A challenging problem lies in the di- 
vergence of the variance of the CR primary flux. As shown re- 
cently by Blasi & Amato (201 1 ), the second moment of the flux 
probability distribution function (PDF) is infinite. Although the 
variance can always be regularized in some way or another, this 
intriguing result threatens the conventional model of CR prop- 
agation, so far extremely successful. That is why we have thor- 
oughly reinvestigated here how the discreteness of the sources 
may modify our vision of Galactic CR propagation. 

To commence, the CR flux of primary cosmic rays obtained 
in the myriad model has a well-defined mean, identical to what 
is derived in the conventional approach by averaging the source 
distribution over space and time. This good point gives credit to 
the continuous model. We have then concentrated on the issue of 
the flux variance which diverges in the myriad model. Although 
the central limit theorem cannot be used in that case, at least in 
its ordinary form, the PDF of the flux is well defined everywhere 
inside the Galactic magnetic halo. We have run Monte Carlo 
simulations with several thousands of different populations of 
point-like sources. For the first time, we have derived the quan- 
tiles of the CR proton flux as a function of proton energy. The 
quantiles turn out to be well-defined and can be used to com- 
pute 68% confidence intervals. In spite of an infinite variance, 
meaningful error bars for the fluxes of primary CR nuclei can be 
defined. These depend on the propagation parameters and tend to 
decrease with the number of sources implied in the signal. That 
is why the uncertainty bands of Fig. [6] widen with the vertical 
extension L of the CR diffusive halo. 

We have so far focused on the flux at the Earth, but the 
same procedure can be applied throughout the Galaxy. It is how- 
ever extremely time consuming and we defer to a subsequent 
publication the analytical derivation of the flux PDF. Such an 
investigation is crucial insofar as secondary species are pro- 
duced by the interactions of primary CR nuclei on the ISM and 
the actual distribution of the latter matters. It would be inter- 
esting to examine the effect of the myriad hypothesis on the 
fluxes of antiprotons or positrons at the Earth as well as on the 
Galactic gamma ray diffuse emission. This emission has been 
so far calculated in the framework of the continuous model. 
An essential prediction is a gamma ray power law spectrum 
which traces the CR proton and helium energy distributions. 
According to the continuous hypothesis, we expect a gamma 
ray spectral index of ~ -2.8 + 0.1 everywhere. Measurements 
by the H.E.S.S. collaboration ( |Aharonian et al.|20 06 ) of the dif- 
fuse emission from the Galactic centre indicate a photon index of 
2.29 + 0.07 + 0.02, significantly below the prediction of the con- 
tinuous model. Fluctuations of the CR proton and helium spectra 
around their mean values have thus already been observed. 

As regards the primary CR fluxes at the Earth, another im- 
portant result is the use of a catalog for the nearby and young 
sources which span the so-called local region. These objects 
have been extracted from the Green and ATNF surveys and are 
based on the latest astronomical observations. The catalog yields 
a contribution <p cat to the total flux <f>, whereas the sources extend- 
ing beyond the local region contribute the complement <f> ext . This 
procedure provides a natural regularization of the flux variance 



since it is based on observations. We have found that below a 
few hundreds of GeV, cat is always negligible with respect to 
0ext, whatever the CR propagation parameters. Since the flux is 
dominated by distant or old sources, the continuous hypothesis 
applies and the use of the conventional CR propagation model 
is fully justified. We have succeeded in reconciling the myriad 
model with the continuous approach. Both should give identical 
results fo the low energy CR fluxes at the Earth. At the TeV scale 
and beyond, the situation becomes more complicated. If the to- 
tal number of objects that source the flux is large with respect 
to the catalog, the contribution (f> ext dominates and the continu- 
ous hypothesis is still valid at high energies. On the contrary, if 
the supernova explosion rate v becomes smaller than the canon- 
ical value of 3 events per century, or if the half- thickness L of 
the DH decreases, the catalog becomes relatively important and 
may even dominate the flux above a few TeV. This is actually 
the case for the black curve of Fig. 11 which corresponds to 



the min model and v = 1 century . The CREAM data points 
lie inside the catalog 68% error band. Notice how well the CR 
proton excess is naturally explained by local and young sources 
which have actually been detected. There is no need for a break 
in the injection spectrum, nor for a peculiar behaviour of the dif- 
fusion coefficient D with energy. It should even be possible to 
find a particuliar set of CR propagation parameters that would 
make the catalog a natural explanation of both PAMELA and 
CREAM data. Since the local region alone is implied, these pa- 
rameters ought not to be the same as for the bulk of the Galactic 
magnetic halo. 



Appendix A: Computation of the Green function 

In order to solve analytically Eq. (p), we need to simplify our 
description of the Milky Way DH and replace it by an infinite 
slab of half-thickness L with a gaseous disk in the middle at 
Z = 0. Radial boundary conditions at r = /? ga i are no longer 
implemented in the propagator. This simplification of the setup 
for CR propagation could be a problem if we were interested in 
the CR densities iff close to the radial boundaries, at a distance 
of 20 kpc from the Galactic center. But our aim is to calculate 
these densities at the Earth, at a galactocentric distance of r Q = 
8.5 kpc, i.e. far from the radial boundaries. Furthermore, even 
though the propagator 8^, is derived within the framework of an 
infinite diffusive slab, integrals on the sources of cosmic rays, 
such as relation (|7]), are still performed up to the radial bound- 
aries at r = Rgd. We have checked that such a procedure does 
not introduce any significant error on the CR fluxes at the Earth. 
In the case of a source term q- dcc that is continuous in space and 
time, this approach yields results very close to those obtained 
with a method based on radial Bessel functions. With this sim- 
plified setup, the propagation of CR species becomes invariant 
under a translation along the horizontal axes x and y. The master 
equation |9]l still needs to be solved along the vertical direction z, 
with the condition that & p vanishes at the boundaries z = +L. The 
construction of the Green function for CR nuclei is inspired from 
the solution to the heat diffusion problem and has been given by 
many authors. The horizontal and vertical dependencies in @ p 
can be factored out by setting 



%(x, f ^x s , fs )=^exp(-^-l 7,,=., 



zs,ts) , 
(A.l) 
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where r = t — t$ and p 2 = (x - x$) 2 + (y - ys) 2 - The vertical 
function 'fp is given by the expansion 

■vc . t , j Vc(\z\-\z s \) \ 
T p (z,t <- zsJs) = exp^ — }xO(t) 



2D 



C' 



(A.2) 



where 0(t) is the Heaviside function. The sum runs over even 
and odd eigenfunctions. The former can be expressed as 



<f„( Z ) = sin{fc„(L~| Z |)} 



(A.3) 



The corresponding even eigenwavevectors k„ are solutions of the 
equation 

- tan(/c„L) = — \- , (A.4) 

where fe^, = V C /2D and /c s = hT p /D are typical wavevectors 
which account for the effects of Galactic convection and proton 
collisions on the ISM. Notice that at high energy, the dimension- 
less parameters k w L and k s L become very small as D increases. 
We then expect convection and spallations to have little effect 
with respect to space diffusion at high energy, as illustrated in 
Fig. |A.1| The odd eigenfunctions are given by 



<(z) = sin{^(L-z)} 



(A.5) 



where the odd eigenwavevectors k' n are such that k' n L = nn. This 
condition ensures that S^(z) vanishes at z — and is an odd func- 
tion of the height z- The decay rates a„ and a' n are respectively 
related to the eigenwavevectors k„ and k' n through 



<*» = Dk 2 n + ^ . 



(A.6) 



a n > = Dk n , + 



V 2 
_£ 

AD 



The eigenfunctions S n (z) and $' n {z) may be understood as or- 
thogonal vectors of a basis over which any vertical function that 
vanishes at z — +L can be expanded. They are normalized in 
such a way that 



<<5 



\ | S n ) = f dz <?/(z) S n {z) = 5 pn C n (A. 



7) 



and 



(&p I O ~ S pn C'„ = S pn L , 

where the normalization factors C„ and C' n are 



C„ = L \ 1 



sinc~(/c„L) 



and C' = L 



(A.8) 



(A.9) 



The dimensionless parameter p is defined by 1 /p — (k w + k s )L. 
It becomes infinite in the absence of Galactic convection and 
collisions on the ISM. In that limit, the eigenwavevectors k„ are 
equal to (« - 1 /2) n/L. 

We still need to address a technical problem related to the 
behaviour of the vertical propagator ( |A.2j > when t is close to t s . 
We shall find that the variance of the CR proton flux is dom- 
inated by local and recent supernova explosions. In that case, 
the period of time r that separates the injection of protons from 
their detection at the Earth becomes small. The expansion \A.2\ 
needs to be pushed a long way up in terms of the number of 



eigenfunctions involved. When r is small, the exponentially de- 
creasing functions exp(-Q', 1 T) and exp(-a^r) are still close to 
unity unless n becomes very large. The numerical convergence 
of y p requires then to sum expression ( A.2 1 over an exceedingly 
large number of terms, hence a potential problem of CPU time. 



Because in that regime, most of the expansion (A.2 1 is provided 



by high order terms, Galactic convection and CR spallations be- 
come negligible with respect to diffusion. For these terms, the 
even eigenwavevectors k n are actually very close to their pure 
diffusion values of (n - 1 /2) n/L, even in the case where p is not 
very large. In this regime of basically pure diffusion, another so- 
lution to CR propagation is provided by the method of electrical 
images. To commence, if the half-thickness L of the slab is made 
infinite, we recover pure diffusion in infinite 3D space. In that 
case, the propagator is well-known, with a vertical contribution 
expressed as 



%(z,t <- z S ,t S ) = fm(z,t 



exp 



ZsJs) 
(z- 



zsY 



4Dt 



(A. 10) 



As already discussed by Baltz & Edsjo (1999J, the vertical 



boundaries of the DH can now be implemented by associating 
to each point-like source lying at height zs the infinite series of 
its multiple images through the planes z = +L. The boundaries 
act as mirrors and give from the source an infinite series of im- 
ages. The n-th image is located at 



2Ln + {-\) n z s 



(All) 



and has a positive or negative contribution depending on whether 
n is an even or odd number. When the diffusion time t = t — ts 
is small, the ID solution (|A.10|) is a quite good approximation. 

iy 



The relevant parameter is actua! 



n = 



ADt 



(A. 12) 



and, in the regime where it is much larger than 1, i.e., for small 
values of r, the propagation is insensitive to the vertical bound- 
aries. When the diffusion parameter rj decreases and the diffusion 
length /Id = V4Dr becomes comparable to the half-thickness L 
of the DH, images of the point-like source need to be taken into 
account in the sum 

%{z,t «- zsJs) = J] rio(z ' r *~ z "' ts) ■ (A ' 13) 

n = -oo 

When r] is much smaller than 1, many terms must be taken into 
account in the previous sum. However, this regime corresponds 



to large values of the diffusion time r for which expansion (A.2 1 
converges very rapidly. 

We have devised two complementary methods to calculate 
the propagator Depending on the value of the diffusion pa- 
rameter 77, we can use the relation ( |A.13[ ) of electrical images 
(77 > 1) or the expansion ( |A.2| i (77 < 1). Since the CR sources and 
the Earth are located inside the Galactic disk, we have set zs and 
z equal to in Fig. |A.l| The vertical part of the Green function 
depends only on the time delay r = t—t$ with "fp(j) = ^(0, t «— 
0, ts) = %(0,t <— 0,0). At low energy, Galactic convection 
and CR collisions on the ISM become important with respect 
to space diffusion. The method of electrical images is reliable 
only for small values of the diffusion time r. Expansion ( A. 1 3] > 
can still be used to calculate y p in the regime where 77 becomes 



9 



G. Bernard et al.: Variance of the Galactic nuclei cosmic ray flux 



large. How large depends on the relative strength of the vari- 
ous CR propagation mechanisms. To get a flavor of the range 
of validity over which electrical images can be used even in the 
case where k s and k w are larger than 1/L, we have borrowed 
the MED set of CR propagation parameters from |Donato et al. 
(2004). This benchmark configuration provides the best fit to the 
B/C measurements (see Table 1). The black short dashed line of 
Fig. |A. 1 [ corresponds to pure diffusion in infinite 3D space. The 
method of electrical images can only be used for pure CR dif- 
fusion and yields the red solid line. Expansion (A. 13 1 has been 
calculated with n sllm = 100 images and is valid down to the small 
value of T] ~ 0.016. Expansion (A. 2 1 has also been pushed up to 



the 100-th term and converges even for 77 as large as 10 . It leads 
to the blue curves which correspond each to a different proton 
kinetic energy T p . Because diffusion takes over the other pro- 
cesses at very high energy, the blue short dashed curve is com- 
pletely superimposed on the red solid line. In that regime, expan- 
sions (A.2i and ( A.13 1 yield the same result. As T p decreases, 



convection and spallations come into play and the blue curves 
depart from their high energy limit. Notice that the 100 GeV con- 
figuration is still fairly close to the high energy case. At 10 GeV, 
the blue dotted-long dashed curve differs noticeably from the red 
solid line. The pure diffusive regime is nevertheless obtained for 
77 > 100. Finally, the 1 GeV blue dotted-short dashed curve is 
significantly shifted towards higher values of 77 with respect to 
the pure diffusive case. In order to reliably calculate the proton 
propagator below ~ 1 GeV, many terms need to be taken into 
account in expansion (A. 2 1. This sum should be used up to large 



values of 77 before the electrical images provide an accurate re- 
sult. 



Appendix B: Probability distribution for the flux 

In this section, we compute the high-0 behaviour of the probabil- 
ity density P((f>) for the flux due to a point source drawn from a 
uniform spatial and temporal distribution. This high-0 tail of the 
distribution is the part leading to the divergence of the variance. 
It is due to the sources which are very close and very young, for 
which the effects of Galactic wind, escape and reacceleration are 
very small. We can neglect these effects, as long as we are only 
interested in the asymptotic behaviour of P(<f>) at high <p. 

The flux at distance r and at time t from a point source emit- 
ting instantly all its particles at r — and t — is given by 



a -P-jAKt 

,3/2 e 



where 



q 



(AnKfl 2 



Case of a 2D distribution of sources (thin disk) 

Let us first consider sources having a given age t . The probability 
density that a source lies at a distance r is given by 



p(r) = 



dP(r) _ 2r 
dr ~ R 2 



where R stands for the radius of the region containing the 
sources. We have 

2r dr 

d(j> 



so that 



p(<f>\t) = 



AKt T 
dP(<f>,t) AKt 1 



d(p R 2 
For an age f, fluxes are in the interval 

a e -R 2 /4Kt <(b< (l 



which can be written as 

AKt 1 la jii 1 ay, a \ 

Where W stands for the window function, being equal to 1 in the 
interval and outside. It is easily checked that 



P(<p\t)d<p 



AKt 

~R2 



In 



1 



The probability distribution for <p is obtained by 

p(0)= J p((j)\t)p{t)dt 



where the upper bound is 

'max — rX 

the lower bound t mm (d>) is a solution of 



,2/3^-2/3 



a 

J/2' 



-R 2 /4Kt m 



If the age distribution pit) is uniform between and T, we have 
p(t) = l/T so that 



P(<P) = 



AK 

R 2 T<p 



J 'rain 



tdt 



p(4>) 



i 2 T(f> 
2Ka 4 ' 3 



R 2 T 

This can be written as 

2Ka A ' 3 



PW 

At high <f>, we have 



R 2 T 



-7/3 



P(<f>)- 



-7/3 _ min vy/ 
rA«4/3 



\\ _ e -R 2 /4»mi»(0)j 



-7/3 



Case of a 3D distribution of sources (thick disk) 

If the sources are distributed in a volume instead of a surface, we 
still have 

rdr 

d<f> = -. 



but 



p(r) = 



2Kt T 

dP(r) _ 3r 2 
dr ~ ~PJ 



dP = 



3r 2 dr _ 6rKt d<f> 
R 3 ~ R 3 ~S 



As 



we have 



= ^-AKt\n(4>t 3 l 2 la) 

_ \2(Ktf 12 V-ln(0f 3/2 /a) 
~ R~ 3 d> 



As before, we can check that the total probability is 1. Let us 
compute 



f3/2 



(3/2 



f 12 ( gf ) 3/2 f V-ln(0f 3 / 2 / a ) 

J*W rf *=— R3— J ^ 



d(f> 
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We set x = (pt 3/2 /a, 



J pith t) ,kh = v "_" - — I - - — — tlx 

The integral is given by 



R 3 



r<t> 

J<t>m 



12(Kt) 3/2 f^" f3/2/a V^ln^ 

,nf 3/2 /a X 



/V- lnx 2 , /2 
dx = - (- In x) 3/z + cte 
x 3 



Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993 
Maurin, D., Donato, R, Taillet, R., & Salati, P. 2001, ApJ, 555, 585 
Nakamura et al., K. 2010, J. Phys. G: Nucl. Part. Phys., 37, 075021 
Norbury, J. W. & Townsend, L. W. 2007, Nuclear Instruments and Methods in 

Physics Research B, 254, 187 
Taillet, R., Salati, P., Maurin, D., Vangioni-Flam, E., & Casse, M. 2004, ApJ, 

609, 173 

Yoon, Y. S., Ahn, H. S., Allison, P. S., et al. 201 1, ApJ, 728, 122 
Yuan, Q., Zhang, B., & Bi, X.-J. 201 1, Phys. Rev. D, 84, 043002 
Yusifov, I. & Kiiciik, I. 2004, A&A, 422, 545 



so that 



P(<p\t)d(p 



8(Kt) 3 ' 2 r , , ^ /9 #max' 3/2 /a 



[i-mx^f 



Finally, as max f 3/2 = a and min f 3/2 = ae~ R2/4Kt , 
f , 1XJ 8(Kt) 3 ' 2 ( R 2 \ 3/2 , 

The probability distribution is obtained by 
P((f>) = J p{0)p{t)dt 



As before, p(t) = l/T, 

P(<P) = f 
ds 

y = tcf> 2 ' 3 /a 2 '\ 



-dt 



\2(Ktf 12 V-ln(0f 3/2 /fl) 
TR 3 ' 

where the bounds f min and f max are the same as before. We define 



pW 



S 



12a 5 ' 3 (Ky) 3 ' 2 yj-\ny 3 l 2 dy 



TR 3 



12a 5 ' 3 K 3 ' 2 
TR 3 ^ 3 



j y 3/2 ^pfoy^dy 



withy max = 1 andy m i n (0) solution of y = exp(-R 2 a(f> 2/3 /6Ky). 
For high values of <p, the lower bound vanishes and the integral 
does not depend on (f>, which yields the final result 

oc <p-^ 3 

Acknowledgements. We thank Prof. Blasi for very useful discussions. This 
work was supported by the Spanish MICINN's Consolider-Ingenio 2010 
Programme under grant CPAN CSD2007-00042. We also thank the support of 
the MICINN under grant FPA2009-08958, the Community of Madrid under 
grant HEPHACOS S2()()9/ESP-1473,and the European Union under the Marie 
Curie-ITN program PITN-GA-2009-237920 



References 

Adriani, O., Barbarino, G. C, Bazilevskaya, G. A., et al. 2011, Science, 332, 69 

Aharonian, F. et al. 2006, Nature, 439, 695 

Baltz, E. A. & Edsjo, J. 1999, Phys. Rev. D, 59, 02351 1 

Blasi, P. & Amato, E. 2011, ArXiv e-prints 

Delahaye, T., Lavalle, J., Lineros, R., Donato, F., & Fornengo, N. 2010, A&A, 
524, A51 

Donato, F., Fornengo, N., Maurin, D., Salati, P., & Taillet, R. 2004, Phys. Rev. D, 
69, 063501 

Donato, R, Maurin, D., Salati, P., et al. 2001, ApJ, 563, 172 

Gradshteyn, I. S., Ryzhik, I. M., Jeffrey, A., & Zwillinger, D. 2007, Table of 

Integrals, Series, and Products, ed. Gradshteyn, I. S., Ryzhik, I. M., Jeffrey, 

A., & Zwillinger, D. 
Green, D. A. 2009, A Catalogue of Galactic Supernova Remnants (2009 March 

version) (Cambridge) 
Higdon, S. J. U., Higdon, J. L., van der Hulst, J. M., & Stacey, G. J. 2003, ApJ, 

592, 161 



11 



G. Bernard et al.: Variance of the Galactic nuclei cosmic ray flux 



\ (4>) 









10 s 
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10 J 
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Fig. 1. Evolution of the mean value and of the standard deviation of the flux with the age cut-off f m j n . The cut-off obtained with the 
prescription f m ; n = R mdx / yjAvD(E) (see text) is indicated with an arrow. This is for the max model, at an energy of 10 TeV. 




Tp (GeV) 



Fig. 2. Proton flux obtained for the min diffusion model (see below) and a source frequency fi = 3 century 1 . The band displays 
the standard deviation cr^ obtained analytically, with a lower cut-off on the ages of the sources f m ; n = 100 years. The blue symbols 
indicate the flux measured by CREAM ( |Yoon et al.|2011| l and PAMELA ( |Adriani et al.|2011) . 




a/R 2a/R 3a/R 4> 

Fig. 3. Probability that the flux tp = a/r is greater than <p, for a unique point source drawn randomly in the disk, within a distance R. 
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Fig. 4. Probability distribution of cp/N for N = 1, 10 and 100 sources. 

P{4>) 105 1 
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Fig. 5. Examples of histograms of fluxes for the "min" cosmic ray propagation model. 
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Tp (GeV) 



Fig. 6. Confidence intervals for the flux from all the sources : the red curves show the 10 % quantiles, the black curves show the 
mean value and the grey band show the 68% confidence interval for the flux, for v = 3 century" 1 . 
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Fig. 7. The bands indicate the theoretical expected number of sources within distance r and for a given age (left), or with ages less 
than t within a given distance (right), for a supernova rate of 3 explosions per century in the Galaxy. The width of the bands gauges 
the uncertainty in the number of local sources due to the shot-noise effect. The curves feature the cumulated number of sources in 
our catalog. This one appears to be complete for t < 5 x 10 4 y and r < 2 kpc. 
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I. 05 





Fig. 9. In each panel, the red lines correspond to the 10% quantiles of the CR proton flux generated by random sources drawn 
in the local region. The mean value is indicated by the black curve. The flux produced by the sources belonging to our catalog is 
represented in blue. The systematic errors associated to the uncertainties on the distances and ages of these objects span the grey 
band (standard deviation). All curves are derived with a supernova explosion rate of v = 3 century -1 . The min, med and max 
propagation benchmark models are considered. 
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TOA proton kinetic energy T p (GcV) TOA proton kinetic energy T p (GeV) 



Fig. 10. Mean flux (black curves) and envelopes representing the standard deviation of the flux, for the min, med and max propaga- 
tion model and for an exponential distribution of sources along z (left) and for a thin disk (right), for v = 3 century 1 . 




Fig. 11. The blue, red and black curves feature the total flux <p computed as the sum of the mean external flux (0 e xt) an d the 
contribution $> cat from the catalog. They respectively correspond to the max, med and min CR propagation benchmark models. The 
bands which extend around the curves have the same meaning as in Fig. [9] They indicate the standard deviation of the flux associated 
to the observational errors on the ages and distances of the SNR of the catalog. 
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Fig.A.l. The effects of Galactic wind and spallations on the vertical propagator ^ are presented in this figure where the MED 
model of CR propagation, which best fits the B/C ratio, has been selected. At very high proton energy, diffusion dominates over the 
other processes and we recover the results obtained by setting k w and k s equal to 0. Because n sum = 100 terms have been taken into 
account in the expansions (A.2i and ( A.13| >, the blue short dashed curve matches exactly the red solid line. As the proton kinetic 
energy T p is decreased down to 1 GeV, Galactic convection and proton interactions on the ISM become more and more important 
and the blue curves departe from their infinite energy limit. 
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